%% Code for simulating paths with optimal contract

%% Monte-Carlo
day_incr = 1;

dt = day_incr/365; %in days
YRS = 50;

T = round(YRS/dt)+1;
n_sim = 10000;

dW = sqrt(dt)*randn(n_sim, T);

W = cumsum(dW,1);

%% Pick prior and hidden states

p0 = 0.9;
n_good = round(n_sim*p0);

G_IDX = zeros(n_sim,1);
B_IDX = zeros(n_sim,1);

hidden_type_on = 1;

if hidden_type_on == 1
    G_IDX(1:n_good) = 1;
    B_IDX = 1- G_IDX;
end


%% INTERPOLANTS

% benchmark
V_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),V_opt(:));

C_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),C_opt(:));
K_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),K_opt(:));
Beta_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),Beta_opt(:));
Eta_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),Eta_opt(:));

ExRet_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),ExRet_opt(:));

mu_v_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),mu_v_opt(:));
mu_y_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),mu_y_opt(:));
vol_y_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),vol_y_opt(:));
vol_phi_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),vol_phi_opt(:));

vol_c_opt_FUN = scatteredInterpolant(Y_mat(:),PHI_mat(:),vol_c_opt(:));


%% GENERATE HISTORIES
y = NaN(n_sim,T+1);
phi = NaN(n_sim, T+1);

y(:,1) = Y(end);
phi(:,1) = p0;

dW_cond = NaN(n_sim,T+1);
dR = NaN(n_sim,T+1);

v = NaN(n_sim, T+1);
v(:,1) = 1;

c = NaN(n_sim,T+1);
k = NaN(n_sim,T+1);
eta = NaN(n_sim,T+1);
b = NaN(n_sim,T+1);
ff = NaN(n_sim,T+1);
mv = NaN(n_sim,T+1);
my = NaN(n_sim,T+1);
sy = NaN(n_sim,T+1);
sphi = NaN(n_sim,T+1);

sc = NaN(n_sim,T+1);

exret = NaN(n_sim,T+1);


% SIMULATION

y_start = y(:,1);
phi_start = phi(:,1);

for t = 1:T
    
    % Controls
    c(:,t) = C_FUN(y_start, phi_start);
    k(:,t) = K_FUN(y_start, phi_start);
    eta(:,t) = Eta_FUN(y_start, phi_start);    
    b(:,t) = Beta_FUN(y_start, phi_start);
    
    sc(:,t) = vol_c_FUN(y_start, phi_start);
    
    % Returns
    exret(:,t) = ExRet_FUN(y_start, phi_start);
    
    % law of motion of state variables
    mv(:,t) = mu_v_FUN(y_start, phi_start);
    my(:,t) = mu_y_FUN(y_start, phi_start);
    sy(:,t) = vol_y_FUN(y_start, phi_start);
    sphi(:,t) = vol_phi_FUN(y_start, phi_start);
    
    %Shocks 
    dW_cond(:,t) = (dW(:,t) + G_IDX.*eta(:,t).*(1 - phi(:,t))*dt + B_IDX.*eta(:,t).*(- phi(:,t))*dt);

    %Realized returns
    dR(:,t) = exret(:,t)*dt + sigma*dW_cond(:,t);
    
    %update state variables  
    y(:,t+1) = y_start + my(:,t)*dt + sy(:,t).*dW_cond(:,t);
    phi(:,t+1) = phi_start + sphi(:,t).*dW_cond(:,t);
    
    v(:,t+1) = v(:,t).*(1 + mv(:,t)*dt + b(:,t).*dW_cond(:,t));
    
    y_start = y(:,t+1);
    phi_start = phi(:,t+1);
    
    %disp(t)
end

elast_c = (sc./c + b)/sigma;

%% ANALYSIS - GOOD

% Mean
y_g_mean = mean(y(logical(G_IDX),:));
phi_g_mean = mean(phi(logical(G_IDX),:));
b_g_mean = mean(b(logical(G_IDX),:));
k_g_mean = mean(k(logical(G_IDX),:));
elast_c_g_mean = mean(elast_c(logical(G_IDX),:));

% Median
y_g_median = median(y(logical(G_IDX),:));
phi_g_median = median(phi(logical(G_IDX),:));
b_g_median = median(b(logical(G_IDX),:));
k_g_median = median(k(logical(G_IDX),:));
elast_c_g_median = median(elast_c(logical(G_IDX),:));

% pc25
y_g_pc25 = prctile(y(logical(G_IDX),:),25);
phi_g_pc25 = prctile(phi(logical(G_IDX),:),25);
b_g_pc25 = prctile(b(logical(G_IDX),:),25);
k_g_pc25 = prctile(k(logical(G_IDX),:),25);
elast_c_g_pc25 = prctile(elast_c(logical(G_IDX),:),25);

% pc75
y_g_pc75 = prctile(y(logical(G_IDX),:),75);
phi_g_pc75 = prctile(phi(logical(G_IDX),:),75);
b_g_pc75 = prctile(b(logical(G_IDX),:),75);
k_g_pc75 = prctile(k(logical(G_IDX),:),75);
elast_c_g_pc75 = prctile(elast_c(logical(G_IDX),:),75);




%% ANALYSIS - Bad

% Mean
y_b_mean = mean(y(logical(B_IDX),:));
phi_b_mean = mean(phi(logical(B_IDX),:));
b_b_mean = mean(b(logical(B_IDX),:));
k_b_mean = mean(k(logical(B_IDX),:));
elast_c_b_mean = mean(elast_c(logical(B_IDX),:));

% Median
y_b_median = median(y(logical(B_IDX),:));
phi_b_median = median(phi(logical(B_IDX),:));
b_b_median = median(b(logical(B_IDX),:));
k_b_median = median(k(logical(B_IDX),:));
elast_c_b_median = median(elast_c(logical(B_IDX),:));

% pc25
y_b_pc25 = prctile(y(logical(B_IDX),:),25);
phi_b_pc25 = prctile(phi(logical(B_IDX),:),25);
b_b_pc25 = prctile(b(logical(B_IDX),:),25);
k_b_pc25 = prctile(k(logical(B_IDX),:),25);
elast_c_b_pc25 = prctile(elast_c(logical(B_IDX),:),25);

% pc75
y_b_pc75 = prctile(y(logical(B_IDX),:),75);
phi_b_pc75 = prctile(phi(logical(B_IDX),:),75);
b_b_pc75 = prctile(b(logical(B_IDX),:),75);
k_b_pc75 = prctile(k(logical(B_IDX),:),75);
elast_c_b_pc75 = prctile(elast_c(logical(B_IDX),:),75);


%% ANALYSIS - ALL

% Mean
y_a_mean = mean(y);
phi_a_mean = mean(phi);
b_a_mean = mean(b);
k_a_mean = mean(k);
elast_c_a_mean = mean(elast_c);

% Median
y_a_median = median(y);
phi_a_median = median(phi);
b_a_median = median(b);
k_a_median = median(k);
elast_c_a_median = median(elast_c);

% pc25
y_a_pc25 = prctile(y,25);
phi_a_pc25 = prctile(phi,25);
b_a_pc25 = prctile(b,25);
k_a_pc25 = prctile(k,25);
elast_c_a_pc25 = prctile(elast_c,25);

% pc75
y_a_pc75 = prctile(y,75);
phi_a_pc75 = prctile(phi,75);
b_a_pc75 = prctile(b,75);
k_a_pc75 = prctile(k,75);
elast_c_a_pc75 = prctile(elast_c,75);


%% PLOTS
paperposition = [0 0 8 4.5];
definition = '-r500';
axis_font = 6;
label_font = 10;
legend_font = 8;

edgecolor = 'k';

linewidth = 1;

color1 = 'b';
color2 = 'r';
color3 = [0.4940 0.1840 0.5560];

marker1 = 'o';
marker2 = '^';
marker3 = '*';

marker_size = 5.5;

marker_spacing = [1:10*365:T+1];

lab1 = 'Average';
lab2 = 'Median';
lab3 = '25th-75th pct';

min_adj = 0;
if p0 == 0.1
    min_adj = 1/10;
end

bytype_pos = 'southwest';
if p0 == 0.1
    bytype_pos = 'northwest';
end


%% Beta - ALL
VAR1 = b_a_mean;
VAR2 = b_a_median;
VAR3 = b_a_pc25;
VAR4 = b_a_pc75;
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3), nanmin(VAR4)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3), nanmax(VAR4)]);
ymin = ymin - (ymax - ymin)*min_adj;
figure;
hold on
l1 = plot([0:dt:YRS+dt], VAR1, 'LineStyle', '-', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l1 = plot([0:dt:YRS+dt], VAR1, 'LineStyle', 'none', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', '-', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', 'none', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot([0:dt:YRS + dt], VAR3, 'LineStyle', '-', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot([0:dt:YRS + dt], VAR3, 'LineStyle', 'none', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l4 = plot([0:dt:YRS + dt], VAR4, 'LineStyle', '-', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l4 = plot([0:dt:YRS + dt], VAR4, 'LineStyle', 'none', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Time (years)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('Incentives ($$\hat\beta$$)','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{lab1,lab2, lab3},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(1).Position = [icons(1).Position(1)*1/4*1.1 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.1 icons(2).Position(2) icons(2).Position(3)];
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.1 icons(3).Position(2) icons(3).Position(3)];
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['incetives_all_ts_p-',num2str(p0*100),'.png'],'-r500');
hold off; clear h; close;



%% Beta - By TYPE
VAR1 = b_g_mean;
VAR2 = b_b_mean;
ymin = min([nanmin(VAR1), nanmin(VAR2)]);
ymax = max([nanmax(VAR1), nanmax(VAR2)]);
figure;
hold on
l1 = plot([0:dt:YRS + dt], VAR1, 'LineStyle', '-', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l1 = plot([0:dt:YRS + dt], VAR1, 'LineStyle', 'none', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', '-', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', 'none', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Time (years)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('Average incentives ($$\hat\beta$$)','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2],{'Skilled','Unskilled'},'Interpreter', 'latex', 'Location',bytype_pos, 'FontSize', legend_font);
icons(1).Position = [icons(1).Position(1)*1/4*1.1 icons(1).Position(2) icons(1).Position(3)];
icons(4).XData = icons(4).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.1 icons(2).Position(2) icons(2).Position(3)];
icons(6).XData = icons(6).XData/4;
legend boxoff
print('-dpng',gcf,['incetives_type_ts_p-',num2str(p0*100),'.png'],'-r500');
hold off; clear h; close;




%% Captal - ALL
VAR1 = k_a_mean;
VAR2 = k_a_median;
VAR3 = k_a_pc25;
VAR4 = k_a_pc75;
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3), nanmin(VAR4)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3), nanmax(VAR4)]);
ymin = ymin - (ymax - ymin)*min_adj;
figure;
hold on
l1 = plot([0:dt:YRS+dt], VAR1, 'LineStyle', '-', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l1 = plot([0:dt:YRS+dt], VAR1, 'LineStyle', 'none', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', '-', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', 'none', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot([0:dt:YRS + dt], VAR3, 'LineStyle', '-', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot([0:dt:YRS + dt], VAR3, 'LineStyle', 'none', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l4 = plot([0:dt:YRS + dt], VAR4, 'LineStyle', '-', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l4 = plot([0:dt:YRS + dt], VAR4, 'LineStyle', 'none', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Time (years)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('Capital ratio ($$k$$)','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{lab1,lab2, lab3},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(1).Position = [icons(1).Position(1)*1/4*1.1 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.1 icons(2).Position(2) icons(2).Position(3)];
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.1 icons(3).Position(2) icons(3).Position(3)];
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['capital_all_ts_p-',num2str(p0*100),'.png'],'-r500');
hold off; clear h; close;


%% Capital - By TYPE
VAR1 = k_g_mean;
VAR2 = k_b_mean;
ymin = min([nanmin(VAR1), nanmin(VAR2)]);
ymax = max([nanmax(VAR1), nanmax(VAR2)]);
figure;
hold on
l1 = plot([0:dt:YRS + dt], VAR1, 'LineStyle', '-', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l1 = plot([0:dt:YRS + dt], VAR1, 'LineStyle', 'none', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', '-', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', 'none', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Time (years)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('Average capital ratio ($$k$$)','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2],{'Skilled','Unskilled'},'Interpreter', 'latex', 'Location',bytype_pos, 'FontSize', legend_font);
icons(1).Position = [icons(1).Position(1)*1/4*1.1 icons(1).Position(2) icons(1).Position(3)];
icons(4).XData = icons(4).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.1 icons(2).Position(2) icons(2).Position(3)];
icons(6).XData = icons(6).XData/4;
legend boxoff
print('-dpng',gcf,['capital_type_ts_p-',num2str(p0*100),'.png'],'-r500');
hold off; clear h; close;




%% PPS - ALL
VAR1 = elast_c_a_mean;
VAR2 = elast_c_a_median;
VAR3 = elast_c_a_pc25;
VAR4 = elast_c_a_pc75;
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3), nanmin(VAR4)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3), nanmax(VAR4)]);
ymin = ymin - (ymax - ymin)*min_adj;
figure;
hold on
l1 = plot([0:dt:YRS+dt], VAR1, 'LineStyle', '-', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l1 = plot([0:dt:YRS+dt], VAR1, 'LineStyle', 'none', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', '-', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', 'none', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot([0:dt:YRS + dt], VAR3, 'LineStyle', '-', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot([0:dt:YRS + dt], VAR3, 'LineStyle', 'none', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l4 = plot([0:dt:YRS + dt], VAR4, 'LineStyle', '-', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l4 = plot([0:dt:YRS + dt], VAR4, 'LineStyle', 'none', 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Time (years)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('Pay-perf. sens. ($$\epsilon_C$$)','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{lab1,lab2, lab3},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(1).Position = [icons(1).Position(1)*1/4*1.1 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.1 icons(2).Position(2) icons(2).Position(3)];
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.1 icons(3).Position(2) icons(3).Position(3)];
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['pps_all_ts_p-',num2str(p0*100),'.png'],'-r500');
hold off; clear h; close;


%% PPS - By TYPE
VAR1 = elast_c_g_mean;
VAR2 = elast_c_b_mean;
ymin = min([nanmin(VAR1), nanmin(VAR2)]);
ymax = max([nanmax(VAR1), nanmax(VAR2)]);
figure;
hold on
l1 = plot([0:dt:YRS + dt], VAR1, 'LineStyle', '-', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l1 = plot([0:dt:YRS + dt], VAR1, 'LineStyle', 'none', 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', '-', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot([0:dt:YRS + dt], VAR2, 'LineStyle', 'none', 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Time (years)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('Avg. pay-perf. sens. ($$\epsilon_C$$)','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2],{'Skilled','Unskilled'},'Interpreter', 'latex', 'Location',bytype_pos, 'FontSize', legend_font);
icons(1).Position = [icons(1).Position(1)*1/4*1.1 icons(1).Position(2) icons(1).Position(3)];
icons(4).XData = icons(4).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.1 icons(2).Position(2) icons(2).Position(3)];
icons(6).XData = icons(6).XData/4;
legend boxoff
print('-dpng',gcf,['pps_type_ts_p-',num2str(p0*100),'.png'],'-r500');
hold off; clear h; close;

